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Psj ■ Abstract 

in : 

We investigate the transition of a quasi-one-dimensional few-boson system from a weakly cor- 
related to a fragmented and finally a fermionized ground state. Our numerically exact analysis, 
based on a multi-configurational method, explores the interplay between different shapes of exter- 
nal and inter-particle forces. Specifically, we demonstrate that the addition of a central barrier to 
an otherwise harmonic trap may supports the system's fragmentation, with a symmetry-induced 

Q-l 

distinction between even and odd atom numbers. Moreover, the impact of inhomogeneous inter- 



^ ■ actions is studied, where the effective coupling strength is spatially modulated. It is laid out how 

the ground state can be displaced in a controlled way depending on the trap and the degree of 
modulation. We present the one- and two-body densities and, beyond that, highlight the role of 
correlations on the basis of the natural occupations. 
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I. INTRODUCTION 



Ever since the first realization of a Bose-Einstein condensate, trapped ultracold atoms 
have been the focal point of an enormous number of research efforts, both from the exper- 
imental and the theoretical side Allowing for an undreamt-of level of control, 
regarding the adjustment of the external as well as the inter-particle forces via electromag- 
netic fields, they may serve as some kind of Rosetta stone for various research areas, ranging 
from solid-state physics to optics. 

A special focus rests on the aspect of a system's dimensionality. Particles restricted to a 
lower-dimensional subspace, such as a wave guide in the one-dimensional case, often unveil 
features that are conspicuously different from isotropic ones. A striking example is the 
scattering theory of an ultracold system whose transverse degrees of freedom are frozen such 
that an effective one-dimensional description becomes possible, as developed by Olshanii 
In that case, the effective interaction strength of the system may be tuned at will from a 
weakly correlated to a strongly repulsive system to an attractive one by merely changing 
the lengthscale of the transverse confinement. A particularly absorbing issue is the so-called 
Tonks-Girardeau gas of impenetrable bosons [(|, which provides an analogy to an ideal gas 



of fermions and has recently been realized experimentally 

Traditionally, the physics of ultracold atoms has been studied extensively on the premise 
of large numbers of atoms N and sufficiently small interaction. This legitimates the use of the 
Gross-Pitaevskii equation which rests on the assumption of a macroscopic wave function 
composed of a single orbital. An extension of this idea to more than one orbital, the so-called 
best mean field has indeed proven to be very efficient in giving a qualitative picture of 
the pathway from condensates to fragmentation, which occurs when the interactions become 
strong enough to deplete the single-orbital 'condensate' [lfll |. 

Still, there are several rationales to consider systems of few atoms, typically N ~ 1 — 
100. For one thing, the interesting situation of a strongly correlated gas is experimentally 
accessible only for small systems. Yet a more fundamental argument is that few-particle 
systems permit a much higher level of control. There is no thermal cloud, as for large 
A, associated with decoherence and energetically dense manifolds of excitations, but a pure 
quantum system. Instead, finite-size effects become relevant, and two-body correlations have 
to be taken into account from the start. Conversely, few-body systems are indeed amenable 
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to ab-initio calculations, making it tempting to analyze their features in detail and without 
resort to uncontrolled approximations. 

While the last statement is true in principle, it is generally highly challenging from a 
computational standpoint. In fact, many attempts to study few- body systems are based 
either on a n a ly tie S0luti0ns for simp ,e m „ d e, SyS ten, UUWU 

or are restricted to 

very few atoms Q- Q| • Moreover, some numerically exact methods have been put 

forward recently. Part of these are actually designed for larger systems but include only few 



orbitals as in double-well traps 
the simple harmonic oscillator 
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HQ 



gneric model systems such as 

2a. 



_while others regard 
or a double well 

The goal of this paper is to investigate the interplay between external forces, on the one 
hand, and the effect of manipulating the inter-particle forces. Our investigation focuses on 
the numerically exact ground state obtained via the Multi-Configuration Time-Dependent 

nnn 

Hartree method [2a, I2jJ- l28| . We consider the example of N = 2, ... ,6 one-dimensional 
repulsive bosons in a double well, whose barrier separating both wells can be adjusted so 
that both a purely harmonic trap as well as large barriers are accessible. We analyze how this 
competes with the effect of an increasing interaction strength, which leads to fragmentation 
and finally fermionization of the ensemble. That interplay is taken one step further by 
considering a setup where the interaction is inhomogeneous, i.e., the inter-particle forces 
depend on the position of a collision, too. This may prove to be a valuable tool on the road 
to extracting single atoms from an ensemble in a controlled way If the interaction 

strength is slightly higher on one side of the trap, the ground state can be shown to be 
displaced to the other side. The nature of this displacement and its dependency on the trap 
configuration as well as the interaction's strength and modulation are studied. 

This article is organized as follows. In Sec. UH the model is introduced and the relevant 
parameter regimes are discussed. Sec. Mil contains a concise introduction to the computa- 
tional method and how it can be applied to our problem. In the subsequent section the 
few-boson system is studied for standard, homogeneous interactions, casting light on the 
passage between the low- and strong-correlation regime, and what effect the trap configu- 
ration has. The one- and two-particle densities are displayed in Sees. IIV Aj / ilVBl while the 
deeper role of fragmentation is highlighted in llVDI The same program is carried out for the 
case of a collisionally inhomogeneous setup in Sec. EJ 
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II. THE MODEL 



In this article we investigate a system of few interacting bosons (JV = 2, . . . , 6) in an 
external trap. These particles, typically atoms with mass M, are taken to be one-dimensional 
(ID) — more precisely, we assume the other two degrees of freedom to be frozen out in a sense 
described below. 

Let the original 3D system be modeled by the Hamiltonian 

H = Y t h i + Y ; V(T i -T j ), 

i i<j 

where h = 2JTP 2 + U(t) is the one-particle (IP) Hamiltonian with a trapping potential U, 
while V is the two-particle interaction potential with certain low-energy scattering param- 
eters ao, r (scattering length and effective range, respectively). It is well-known that for 
sufficiently low momenta (Axo <C 1) the details of the interaction become irrelevant. More 
precisely, if any other interaction potential is taken, the energy as well as the asymptotic 
wave function will remain unchanged as long as a agrees (and r , in the next order). In 
particular, it is usually convenient to model V by an effective point interaction, the so-called 
regularized pseudo-potential j^lj 

V(r) = ^5{v)d r r. 

A. Effective ID Hamiltonian 

We are now interested in the quasi-lD case, where the trap is supposed to be anisotropic 
such that there is a 'transversal' direction (_L) with a characteristic length a± much smaller 
than that of the longitudinal direction, an. In other words, the transverse energy gap be 
sufficiently large compared to the accessible energy of the system, so that (||) may be regarded 
as virtually the only degree of freedom. In this case, one may integrate out the 'frozen' 
transversal subsystem so as to obtain an effective ID interaction \^ 

V(x) = g m 5(x), with g m = ^-(l- C^-) (C = 1.4603 . . . ). 

Ma ± \ a_|_ J 

(These results were derived on the premise of a harmonic transverse trap potential, that 



is, a± = l/y/Mcj±.) It is this effective interaction that shall serve as the base for our 
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investigations, yet with two qualifications. Firstly, the homogeneity of V will be aban- 
doned later in Sec. El where that modification is discussed in detail. Secondly, though this 
delta-type potential is convenient and yields an immediate connection to the experimen- 
tally relevant parameters, it is notoriously intractable in ab-initio computations like ours. 
This is a fundamental fact: By construction, the pseudopotential imposes the condition 
^'(0 + ) — i/j'(Q~) = 2(yf 1D ^(0) on the derivative of the relative coordinate — i.e., whenever two 
particles meet, the wave function behaves like e~ K \ x \. This non-smoothness is of course un- 
physical and solely serves to impose the correct asymptotics on the wave function. In an 
exact calculation, when the problem is approximated by C°°-functions, this leads to conver- 
gence problems. It would be much less artificial to use a more realistic interaction with a 
non-zero effective range as a remedy. We thus opt to mollify the delta function, and instead 
use V(x) = gm5 a {x) with the normalized Gaussian 



model calculations, we ascertained that for a <C 1/|<7id| the results are actually quite close 
to the limit a — > 0. On the other hand, the width a should not be too small so as to 
accurately sample the Gaussian. As a trade-off, we choose a fixed value cr/ay = .05. In 
principle, computations for more than one width could be done in order to extrapolate to 
zero; however, we will always keep it fixed. 

B. Scaling 

For reasons of universality as well as computational aspects, we will work with a Hamil- 
tonian rescaled to the lengthscale of the ID-longitudinal system, ay. More specifically, we 
carry out a global coordinate transform Q' := Q/au, with Q = (x\, . . . , Xjv) T , which leads to 



Here u)\\ = 1/Ma| defines the energy scale, and U'(x') := U(x = x'a\\)/uj\\ etc. denotes the 
rescaled potential deprived of any dimensionful parameters. H' naturally lends itself as a 
convenient working Hamiltonian, and we will skip any primes in the following sections. 




which tends to S as a 



in the distribution sense. In both analytic and numerical 
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As an illustration, the ID point interaction reduces to 




) 



i 



(1) 



The only relevant parameter of the interaction is thus the scaled interaction strength, which 
in turn requires only the knowledge of the (scaled) scattering length a' = a /a\\ and the 
transverse dimension a' ± = a±/a\\. 

C. Parameter regimes 

As mentioned above, two parameters enter our Hamiltonian: a' = ao/ay and a' ± = a±/a\\. 
Both of course depend on 

• the ID length scale a\\ = 1/ yjMu\\ (due to scaling) 

• the scattering length a < an of the atomic species considered (about 100 a.u. for 
alkalis; only positive values are considered here). 

• the transversal lengthscale a± <C an . Of course a± > a is required unless the validity 
of the 'bare' pseudopotential is put into question. We put a± = O.lay for simplicity. 

According to (JTJ> , g\n does not depend linearly on a , but rather tends to +oo as a — > a±/C 
from below. In other words, the system becomes strongly correlated when the scattering 
length approaches the transverse-confinement scale, no matter if the 3D system was strongly 
interacting to begin with. We restrict our attention to g = gxo > 0. Table |U illustrates the 
range of values of a' Q for different (longitudinal) trap frequencies u>\\, and what g' m they 
correspond to for Na/Rb (at fixed a' ± = .1). 



III. COMPUTATIONAL METHOD 

Our goal is to investigate the ground state of the system introduced in Sec. HU for all 
relevant interaction strengths in a numerically exact fashion. In other words, our approach 
is not to approximate the problem by resorting to two-mode or mean-field descriptions, but 
rather to approach the solution in a controllable way. It has to be stressed that this is a highly 
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Table I: Values of the scaled coupling strength g' 1D for Sodium and Rubidium for different trap 
frequencies un/2-n and a' ± = .1. 

challenging and time-consuming task, and only few such studies on ultracold atoms exist 

nnn 

even for model systems I2fj, |21j, 122 ■ Our approach relies on the Multi-Configuration Time- 

nnn 

Dependent Hartree (mctdh) method [26, 28, l22fl, primarily a wave-packet dynamics code 
known for its outstanding efficiency in high-dimensional applications. To be self-contained, 
a concise introduction to this tool — and how it can be adapted to our purposes — is presented 
in this section. 

The underlying idea of MCTDH is to solve the time-dependent Schrodinger equation 



(2) 



*(Q,0) = *„(<?) 

as an initial- value problem by expansion in terms of direct (or Hartree) products $j: 

n\ n f f 

*(Q,t) = j2 A AmAQ,t) = J2 ■ ■ ■ E A h ... jf (t)i[^(x K ,t), (3) 

J ii=i i/=i k=i 

using a convenient multi-index notation for the configurations, J = (ji---jf), where / 

denotes the number of degrees of freedom and Q = (xi, . . . , xf) T ■ The (unknown) single- 
In) 

particle functions (pjj are in turn represented in a fixed, primitive basis implemented on 
a grid. For indistinguishable particles as in our case, the single-particle functions for each 
degree k — 1, . . . , N are of course identical in both type and number (<fj K , with j K < n). 

Note that in the above expansion, not only the coefficients Aj are time-dependent, but 
so are the Hartree products $j. Using the Dirac-Frenkel variational principle, one can de- 
rive equations of motion for both Aj, $j Q|. Integrating this differential-equation system 
allows one to obtain the time evolution of the system via (JH]). Let us emphasize that the 
conceptual complication above offers an enormous advantage: the basis {$j(-,t)} is varia- 
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tionally optimal at each time t. Thus it can be kept fairly small, rendering the procedure 
very efficient. 

It goes without saying that the basis set above is not inherently permutation symmetric, 
as would be an obvious demand when dealing with bosons. However, the symmetry can be 
enforced by symmetrizing the coefficients Aj, even though this turns out to be unnecessary 
as far as the ground state is concerned, which is automatically bosonic jjj^. 

The Heidelberg mctdh package (2J|, which we use, incorporates a significant extension 
to the basic concept outlined so far, the so-called relaxation method j^fj. mctdh provides 
a way to not only propagate a wave packet, but also to obtain the lowest eigenstates of the 
system. The underlying idea is to propagate some wave function \l>o by the non-unitary e~ Hr 
(propagation in imaginary time.) As r — > oo, this automatically damps out any contribution 
but that stemming from the true ground state |0), 

e- Sr * = TT i e- E < r \J)(J\*o). 
j 

In practice, one relies on a more sophisticated scheme termed improved relaxation. Here 
(ty\H — E\ty) is minimized with respect to both the coefficients Aj and the configurations 
$j. The equations of motion thus obtained are then solved iteratively by first solving 
for Aj(t) (by diagonalization of ((&j\H\§ K )) with fixed $j) and then propagating <3>j in 
imaginary time over a short period. The cycle will then be repeated. 

As it stands, the effort of this method scales exponentially with the number of degrees 
of freedom, n N . Just as an illustration, using 15 orbitals and A = 5 requires 7.6 • 10 5 
configurations J. This restricts our analysis in the current setup to about A = O(10), 
depending on how decisive correlation effects are. If these are indeed essential, then it will 
turn out later that at least n = A orbitals are needed for qualitative convergence alone, while 
the true behavior may necessitate about 15. By contrast, the dependence on the primitive 
basis, and thus on the grid points, is not as severe. In our case, the grid spacing should of 
course be small enough to sample the interaction potential, and we consider a basis set of 
75 harmonic-oscillator functions. 
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U(x) 




-3-2-1 12 3 



Figure 1: Sketch of the model potential U(x) = \x 2 + h5 w (x), consisting of a harmonic trap plus 
a normalized Gaussian of width w = 0.5 and barrier strengths h = 2,5, 10. 

IV. BOSONS IN A DOUBLE WELL 



In this as well as in the following section, we consider the ground-state properties of 
bosons in a double-well trap modeled by 

U(x) = ^x 2 + hS w (x). 

This potential is a superposition of a harmonic oscillator (HO), which it equals asymptoti- 
cally, and a central barrier which splits the trap into two fragments (Fig. QJ. The barrier is 
shaped as a normalized Gaussian 5 W of width w and 'barrier strength' h. 

As w — > 0, the effect of the barrier reduces to that of a mere boundary condition (since 
$w S)i an d the corresponding one-particle problem can be solved analytically |l2l . Isfil ]. 
Although this soluble borderline case presents a neat calibration, the exact width w does 
not play a decisive role, as long as it is larger than the grid spacing and w < 1 so as to 
confine the barrier's effect to the central region. We choose w = 0.5 as a trade-off. 

For h — 0, the case of interacting bosons in a harmonic trap is reproduced. In Sec. IIV Al 
we witness the transition from a simple, weakly interacting condensate (g — > 0) to fragmen- 
tation and finally the Tonks- Girardeau limit (g —>■ oo). As h —>■ oo, the energy barrier will 
greatly exceed the energy available to the atoms, and we end up with two isolated wells. 
Higher g then affect only the fragmentation within each of these wells. In between, there is 
an interesting interplay between the 'static' barrier (h) and 'dynamical barriers' in the form 
of inter-particle forces (g). We study this intermediate regime on the examples of h G {2, 5} 
in Sec. ITVTfl 
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A. The reference case: h = 



In the absence of a central barrier, we simply deal with a harmonic trap, which in the 
case of iV = 2 and point interactions (a — > 0) constitutes an exactly solvable problem 
Starting with all interactions turned off (g = 0), the ground-state solution is simply the 
uncorrelated product of the HO-ground state, \I/ = 4>® N . As long as g is small enough, this 
behavior remains qualitatively unaffected by interactions, which only alter the shape of the 
single-particle functions via a mean field (the Gross-Pitaevskii regime). Increasing g has 
the effect of continually depleting the interaction regions {xi = Xj}; the repulsion forces 
the atoms to isolate each other, a process referred to as fragmentation (cf. Sec. IIV Dl for a 
rigorous account of this). Driven to extremes as g — > oo, this fragmentation saturates in a so- 
called fermionized state. In that case the Bose- Fermi mapping [6] asserts that letting g — > oo 
emulates the effects of the Pauli exclusion principle, and the bosons have accomplished to 
minimize their density overlap. 

These qualitative considerations materialize in the reduced densities of the ground-state. 
In the 2-particle density ^2(^1,^2) — the diagonal kernel of the reduced density operator 

p 2 := tr 3 .jv|#)(tf|, 

which yields the probability density of simultaneously finding any two particles at X\ and 
xi — the correlation diagonal {x\ = x 2 } forms an ever deeper dip (the 'correlation hole'). 
This is illustrated in Fig. [2] for g = 0.2 and g = 194. These represent the borderline cases of 
a weakly interacting system and the strong-correlation regime, respectively. In the former 
case, the density is simply Gaussian-like. For large enough repulsion, in turn, the density 
resembles more and more the checkerboard pattern produced by a (polarized) fermionic 
state. 

The 1-particle density pi(x) offers another tool to visualize the fermionization process. 
Fig. OH gives an impression of how the profile changes from a harmonic one {g <C 1), to one 
flattened due to repulsion for mediate g, and finally to the Tonks-Girardeau limit (g > 15), 
where N humps emerge, mimicking the fermionic behavior. 
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-3-2-10123 -3-2-10123 

X 1 X 1 

(a) (b) 
Figure 2: 2-particle density p2{x\,X2) for a harmonic trap (N = 6) for (a) g = 0.2 and (b) g = 194. 




Figure 3: 1-particle density pi(x) for a harmonic trap (N = 5) for different interactions gir>. 
Note how the profile changes from a weakly interacting one (g = 0.2) to a flattened one due to 
fragmentation, and finally to a fermionized profile featuring N humps (g > 15). 

B. Central barrier h > 

We now introduce a central barrier h5 w (x), so as to turn the harmonic trap into a double 
well. Upon increasing g, there are now two competing effects, exemplified on the 2-particle 
density (Fig. BJ: For small enough g, the barrier h dispels the atoms from the center x = 0. 
This state is uncorrelated: they are localized in both wells regardless of whether there 
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X 1 X 1 

(a) (b) 

Figure 4: 2-particle density ^2(^1,^2) for a double-well trap with barrier strength h = 5 (N = 6). 
(a) The weakly interacting limit g = 0.2 (b) The complementary case of fermionization, g = 194. 

are already any other ones. Indeed, in Fig. |4ja) the diagonal {x\ = £2} shows only a 
tiny depletion due to repulsion. In the complementary case g — > 00, the barrier will be 
almost outweighed by the inter-particle repulsion. In this Tonks-Girardeau limit, the atoms 
distribute so as to minimize their density overlap, at the price of an additional potential 
energy in the barrier region. Hence the two-particle density in Fig. EJb) differs from the 
harmonic counterpart only in the stripes along X1/2 = 0, which indicate suppression in the 
central-barrier region. The question as to what happens in the intermediate region requires 
a distinction between even and odd particle numbers; it will be the focus of the ensuing 
paragraphs. 



Even N : assisted fragmentation 

For an even number of atoms — in what follows iV G {2,4,6} — they initially (at g — 0) 
populate each well with N/2 atoms. Upon raising the interaction strength, they seek to 
isolate each other, which can best be done by having a fragmentation within each well, 
interfering only little with the central barrier. This situation is depicted in Fig. El for four 
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(a) (b) 

Figure 5: 1-particle density in a double well with h = 5: Even vs. odd number of atoms, (a) For 
N = 4 atoms, the fragmentation process is in part supported by the central barrier. By contrast, 
(b) indicates in what sense fermionization is suppressed for an odd number (N = 5), since by 
symmetry one particle should eventually reside in the center. 

atoms. The plots for g = 4.7, 15 reflect the tipping from a coherent state to fragmentation. 
For high enough g, the accessible energy approaches that of the barrier, Nhj '\^2irw, so 
that tunneling becomes more and more dominant. Beyond that point, the density profile 
is expected to resemble that of a purely harmonic trap with slightly suppressed amplitude 
in the barrier region. In particular, for any finite h, the Tonks-Girardeau limit will be of a 
generic form in that it exhibits N density maxima. In conclusion, one might say that the 
fragmentation here is assisted insofar as the central barrier helps isolate the ensemble, a 
statement put more precisely in llVDl 

Odd N: delayed fragmentation 

In the case of an odd number (N = 3,5), the situation differs. At g = 0, the atoms 
are again coherently distributed over the two wells. As we strengthen the interactions, they 
try to enlarge their distance — but by symmetry, this can now only be done by placing one 
particle at x = 0, which in turn is impeded by the barrier (cf. Fig. lEb). In other words, 
the system will have to pay the added interaction energy and distribute the extra particle 
over the two wells, until the former one becomes high enough to afford the place in the 
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central-barrier zone. In that sense, the fragmentation is attenuated by the double- well trap. 

C. Ground-state energy 

One central aspect of our analyses is the ground-state energy. Fig. M depicts a typical 
evolution, E(g), as a function of the coupling strength. As suggested above, the g — > oo 
state is isomorphic to a non-interacting fermionic state. In this light, also the ground-state 
energy E(g) (Fig. EJ) may be interpreted as connecting the free bosonic (g = 0) and the free 
fermionic value, corresponding to the saturation as g — > oo. 

The effect of the interaction at g = can be measured by the slope 

dE MiV-l), Ip . . a ^ N(N - 1) f. , 
^"(0) = -{OOl^^i -x 2 )|00) i J |0o(x)| 4 rfx, 

given by the density overlap of two atoms in the non-interacting ground state. The cen- 
tered harmonic-oscillator orbital 0ho by construction has a low curvature (i.e., kinetic 
energy), thus producing a rather high density overlap. It is thus more susceptible to 
the onset of interactions. By contrast, the presence of a central potential-energy bar- 
rier (h — > oo) evokes an orbital 0dw delocalized in both wells. Its density overlap in 
turn will be smaller, which can be seen schematically by assuming for a moment that 
<t*Dw( x ) ~ ^75 [(puo( x ~ x o) + <Pro( x + x o)] is built from a HO orbital centered in both min- 
ima ±x . Neglecting the density overlap between the right- and left-hand contributions, 
/ |0dw| 4 — I / |0ho| 4 5 suggesting that in a double well, the atoms will feel a lesser effect 
when interactions are turned on. This can be seen in Fig. H3 

D. The role of fragmentation 

We have so far relied on a more intuitive notion of fragmentation. It is natural to ask 
whether some of our previous assertions can be put more quantitatively. This we seek to do 
in the present subsection, not only to highlight some deeper concepts, but also to show why 
a numerically exact approach is so vital. 

It has been argued that ever stronger interactions introduce correlations to the system, 
viz., the one-particle Hamiltonian will be no longer dominant, but the influence of the 
two-particle operator V takes over. The latter one imprints explicit two-body correlation 
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h=2 
h=0 (HO) 



60 80 100 120 140 160 
coupling strength g 1D 



Figure 6: Energy E(g) for the case N = 3. Note the slightly different effect of the interaction, 
measured by the slope at g = 0, for different barrier strengths h = (harmonic trap) and h = 2,5. 
The saturation as g — > oo corresponds about to a fermionized state. 

terms on the eigenvector, as illustrated, for instance, on the example of two atoms in a 



37T | starting from the 



trap | 111 LL2j or the fermionization limit g — > oo solvable for any N 
Bose-Fermi mapping, 

oc e- |Q|2/2 Y[ \xi-Xj\- 

l<i<j<N 

Technically speaking, it is therefore far from clear how long \l/ can be well approximated 
by states composed of one-particle functions, let alone a single such configuration as in 
mean-field approaches (see [10] and references therein). 

Being exact, our method offers a handle on these convergence questions, based on a 
criterion put forward by Penrose Consider the spectral decomposition of the one- 

particle density matrix 

Px = tr 2 .jf\*){y\ = '52Tia\<f>*)(<f>a\, (4) 

a 

where n a G [0, 1] is said to be the population of the natural orbital (f) a . Obviously, if all 
n' a = n a N G N (J2 a n 'a = N), then the density may be mapped to the (non-interacting) 
number state jng,^, . . . ) based on the natural one-particle basis; for non-integer values it 
extends that concept. In particular, the highest such occupation, n , may serve as a measure 
of non-fragmentation: for n = 1, a simple condensate is recovered. This is the well-known 
borderline case of the Gross-Pitaevskii eq.: as g — * 0, p\ — > |0o)(0o| 13- The complemen- 
tary fermionization limit (g — > oo) has been investigated semi-analytically drawing on the 
Fermi/Bose mapping, yielding n ~ iV~- 41 for a harmonic trap 
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Figure 7: Maximum natural population no (g) for (a) a harmonic trap (h = 0) and (b) a double 
well (h = 5). For even N, the fragmentation is enhanced by the higher barrier, whereas for odd N 
it is slowed. 

In between those extremes, a thorough many-body treatment is indispensable. This is 
illustrated in Fig.0 where n (g) is plotted. Observe that, for one thing, the high-correlation 
value n (g — ► oo) is much too low for N = 6 compared with the above value n ~ .48, a 
numerical effect explained below. Moreover, the lines for N = 4 and N = 5 in the double- 
well trap (Fig. 0>) are partly reversed as compared to the harmonic trap, while the N = 3 
value of n is shifted upward. This reflects the suppression of fragmentation for odd atom 
numbers with respect to even N. 

The role played by correlations may also be highlighted by the number of single-particle 
functions needed in a calculation, n (cf. IIIII) . It also determines how many terms in the 
spectral expansion of pi (Eq. are included. Figures |HJa-c) show the convergence of the 
ground-state energy of N = 6 bosons as n is varied. For rather low g = .406, the convergence 
for the harmonic trap is fairly smooth, and the variations are altogether relatively small. 
This asserts that Gross-Pitaevskii works qualitatively well, though it may already take many 
orbitals to achieve a high accuracy. For a double well with h = 5, the convergence is much 
more abrupt. Adding just another orbital, n = 2, lowers the energy drastically in comparison 
with any further refinements. This is intelligible, regarding the fact that two-mode models 
based essentially on the anti-/symmetric orbitals of a double well are widely used in that 
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Figure 8: Convergence of the energy E{n) in the number of single-particle functions n for six bosons. 
For weak interaction g = .406, the difference between the harmonic (a) and the split trap (b) is 
conspicuous, though in both cases Gross-Pitaevskii (n = 1) is not too far off. As opposed to that, 
the onset of fermionization (g = 15) wipes out that difference up to a residual energy offset due to 
the central barrier (c). Here n = N orbitals are needed to be at least qualitatively correct. 

context. 

Fig. [He) in turn casts a light on the strong-correlation regime, g = 15 (in fact its ground 
state for h = 5 is already atop the central energy barrier.) Here the lines E{n) for h — 0, 5 
are almost parallel, shifted only by the energy offset of about h ~ J hS w (x)dx. Clearly 
simple mean-field theory is off by more than a factor of 2, as it drastically overestimates 
the interaction energy and fails to reproduce the characteristic Tonks-Girardeau profile with 
iV humps. Just when n = N = 6, the convergence suddenly settles, and the qualitatively 
correct behavior can be observed. Interestingly enough, adding another orbital, n = 7, has 
virtually no effect on that energy scale: the first N natural occupations all turn out to be 
n = 1/N, corresponding to a number state |lo, . . . , Ijv-i) m the natural-orbital basis. This 
is of course quantitatively incorrect in light of the expected behavior n ~ N~ 0A1 , which 
really can be recovered if n > 15 single-particle functions are included. But apparently that 
ground state seems to be a somewhat stable intermediary solution in the fermionization 
limit if the subspace span{$ j | _L < n} is taken too narrow. This provides an insightful link 
to multi-orbital approaches joL llfl|. 

Let us remark that the criterion employed above is conventional and amenable, but not 
imperative. Indeed, the two-particle density pi by nature reveals correlation effects even 
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more clearly. 



V. INHOMOGENEOUS INTERACTIONS 



We have so far relied on the assumption of homogeneous two-particle forces. These are 
invariant under global translations and thus depend on Xi — Xj alone. While this premise is 
most natural from a fundamental point of view, we should keep in mind that our description 
is not a fully microscopic one, even if we ignore the internal structure of the underlying atoms. 
Rather, it is an effective model stripped not only of the transverse degrees of freedom, but 
of course also of the electromagnetic fields that manipulate both external and inter-particle 
forces. 

With this in mind, it appears legitimate to conceive situations where the strength of the 
interactions depends in addition on the position where the collision takes place, as was done 
in a mean-field framework in Ref. 2lJ (see also citations therein). This may be induced 
by means of a Feshbach resonance, tuning ao(B) by adding a spatial dependence to the 
magnetic field. In our one-dimensional setting, it seems even more convenient to exploit the 
parametric dependence on the transverse subsystem, and modify a± locally so as to imprint 
a spatial dependence on gi^. 

Without reference to the specific realization and its concrete experimental constraints, we 
perform a case study where a generic model for the inhomogeneity is assumed to begin with. 
This model will be presented in Sec. IV Al The interplay of that dynamical inhomogeneity 
with the external forces will be studied for a harmonic (|V BJ) and a double- well trap (|V CJl . 



A. Model interaction 



Whereas modeling a position-dependent interaction in a mean-field description (as in 
is straightforward, since one only has an effective one-particle problem, one faces a 
conceptual problem when using a many-body framework. In general, the coupling would 
depend on both participants Xi,Xj, which is technically possible if somewhat awkward. For 
it to make sense intuitively, we require that its modulation lengthscale be much larger than 
the 'radius' of collision, a. 

With this is mind, it is natural to model our interaction in terms of the respective relative 
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coordinate r := X( — Xj (for fixed and — in order to keep V formally symmetric — the 



center of mass 2R := X{ + Xf 



V(r,R)=g(R)5 a (r). 
The demand that g change slowly over the support of V can be cast as 



> a 



(5) 



in the supremum norm. 

There are various possibilities just what scenario should be examined, be it some kind of 
collision-enhanced tunneling or dynamical self-trapping 39!]. We concentrate on a specific 
model where g is essentially imbalanced between the right- and left-hand sides of the trap 
(Fig. E) : 



g{R) = go 



1 + a tanh 



R 



This signifies that for \R\ 3> L, the coupling takes on the asymptotic values 

g± = lim g(R) = g (l ± a), 

\R\— >±oo 

while it changes on a scale of L near the trap's center about go. The parameter a regulates 
both the relative difference between the asymptotic strengths and their ratio: 



As 

9+ 
9- 



\g± - go\ = go®, 

l + a 



a 



Eq. ((HJ) can be met if L 3> era, which is effortlessly fulfilled if we choose L = 1 for convenience. 



B. The reference case: h = 



Generally speaking, the ground state of atoms immersed in a harmonic trap will be 
centered at the trap's bottom, assuming that we start with a weakly interacting ensemble. 
Hence the modulation of the coupling strength g beyond the center will pass them largely 
unnoticed. It is only for strong enough repulsive interaction — where fragmentation sets in — 
that the density profile will start to split and shift partly outward, thus experiencing an 
asymmetry. 

This picture is supported by our calculations, as demonstrated in Figure [101 for N = 5 
atoms. For low enough go = .4, measuring the average interaction strength, the harmonic 
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Figure 9: Our model of the position-dependent coupling g(R)/go = 1 + atanh (j;). The relative 
modulation, here a = 0.5, determines the asymptotic difference from the average value go, while 
the modulation length L = 1 shall remain fixed. 




X 



Figure 10: 1-particle densities pi(x) for a harmonic trap (N = 5) in the case of inhomogeneous 
interactions, here a = .5. The profile features an imbalance for smaller interactions go, where the 
wave packet is centered too much to sample the modulation of g(R). When fragmentation sets in, 
the profile splits and the asymmetry becomes more distinct. In the fermionization limit, the energy 
costs of an imbalance become too large to keep it up. 

profile is barely altered from the homogeneous case a = 0. An imbalance is noticed for 
medium g = 4.7: the atoms are now able to sample the modulation of the coupling strength 
and find it more inexpensive to locate in the less repulsive zone x < (#-). However, this 
effect ceases as the repulsion becomes larger (g > 15). This may be interpreted as follows: 
the energetical costs for concentrating several particles near one spot are soaring, and this 
in total eventually outweighs the relative energy savings reached by an imbalance. 

A look into the two-body correlations P2(xi, x 2 ) in Fig.[TT]helps us clarify what happens. 



21 



-3-2-10123 -3-2-10123 

X 1 X 1 



(a) (b) 

Figure 11: 2-particle density for a harmonic trap in the presence of inhomogeneous interactions 
(N = 5). (a) For g = 0.4, the packet is localized about the center, thus widely ignoring the 
modulation, (b) When fragmentation sets in (g = 15), it starts to delocalize and consequently 
shifts to R < 0. 

Along the correlation diagonal {x\ = x 2 }, R = %i/2 holds. It is here that the modulation 
can have an impact, whereas for X\ = —x 2 , g(R) = g(0) = go as usual. Clearly the density 
on the diagonal {x\ = x 2 } must be spread enough for the modulation to become effective. 
This is not the case for small interactions. Indeed, for g = 0.4, the packet is localized 
about the center, thus widely ignoring the modulation. Yet for medium g (Fig. HTb) . the 
repulsion-driven broadening has become distinct enough so that the ground states exhibits 
some left-right asymmetry. For strong fragmentation, the correlation diagonal will in turn 
be fully depleted, so obviously the atoms can no longer realize the modulation, and the very 
premise of an inhomogeneity-based displacement has become obsolete. 

The above findings are nicely wrapped up in Fig. [TJl showing graphs of (x) = tr(pix) 
as a function of g for iV = 5. For a = 0, and of course for g = 0, no modulation exists 
and, by symmetry, (x) = 0. Notably, the same goes for g — ► 00, when the correlation 
diagonal is depleted as delineated above, even though the displacement will vanish only very 
slowly. There is a trade-off in between for which (x) becomes extremal. The value where 
this occurs, g^a), depends only weakly on the relative modulation a — despite the fact that 
the maximum ground-state displacement — (x)* will of course increase monotonically with 
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Figure 12: The ground-state displacement —(x) as a function of the average interaction go (N = 5). 
Its universal behavior is an increase up to a maximum value followed by a slow decay. The increase 
at go = is strongly enhanced in the presence of a barrier h > 0, while for the purely harmonic 
trap [h = 0), it is rather slow. Of course the maximum itself is much more pronounced for higher 
modulations a, while being absent in the homogeneous case a = 0. 

a. 

C. Central barrier h > 

In the presence of a sufficiently strong barrier, the situation is a different one. To begin 
with (g = 0), the atoms are not centered as before but rather coherently distributed over the 
two wells. Hence, upon switching on the inhomogeneous interaction, they can immediately 
feel the full impact of its modulation on both sides. For finite barrier strength h, they can 
then re-distribute so as to find a compromise between minimum repulsion and potential 
energy. Even though this mechanism is universal for any particle number, we will lay it out 
for both even and odd N so as to keep a link to the reference case a = 0. 

Even N 

The above process is illustrated in Fig. H31 which evidences an immediate shift from 
the right well to the left one, where the repulsion is weaker. This still corresponds to the 
Gross-Pitaevskii regime of a single dominant orbital: there is no correlation hole; in fact 
the probability density of finding both particles in the left well, p2(—xo, —Xo), may even be 



23 




-3-2-10123 -3-2-10123 

X 1 X 1 

(a) (b) 

Figure 13: 2-particle density for N = 4 bosons in a double- well trap (h = 5) and with inhomogeneous 
interactions (a = .5). (a) Already for g = 0.2, the probability of finding any two atoms in the left 
well is significantly enhanced, (b) At the onset of fragmentation (g = 15), the diagonal {x\ = x 2 } 
is being depleted. 

larger than that for separation, p 2 (±xo, =F£o)- As the interaction passes a critical strength, 
fragmentation sets in, somewhat more pronounced on the right-hand side (Fig. IT31b). Note 
how the diagonal {x\ = x 2 } is being emptied, signifying the incipient destruction of the 
imbalance. 

This reflects in the 1-particle density displayed in Fig. HU The density is almost 'in- 
stantaneously' shuffled from the right to the left. In the curve for g = 4.7, it becomes 
apparent that the fragmentation essentially kicks in separately for both wells, where only 
the right well exhibits the typical repulsion-induced split-up. As asserted already for the 
harmonic reference case (h — 0), the modulation becomes marginal relative to the overall 
fermionization process. This may also be discerned here: for an even number of bosons, the 
asymmetric effect fades, and as in Sec. IIVB| the fragmentation is assisted by the central 
barrier insofar as it now enters separately in the two wells. 
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(a) Four atoms 



(b) Five atoms 



Figure 14: 1-particle density for a double well (h = 5) and modulated coupling strength (a = .5) : 
Even vs. odd number of atoms. In both cases, fragmentation essentially sets in separately in each 
well (e.g., see go = 4.7). However, for N = 4 (a), the fermionization process is again supported by 
the central barrier, while (b) suggests that for N = 5 the extra particle can now be temporarily 
accommodated by the left well. 

OddN 

The behavior here is wholly analogous to that evidenced in Fig. [T3J On the correlation 
diagonal, the density again experiences an asymmetry even for tiny g > 0, where the mean- 
field behavior is still virulent. For a strong enough modulation, say a = .5, coincidence 
of two atoms in the left well (xi/2 = — xo) is even enhanced with respect to separation 
{x\ = —X2), while coincidence in the right well is extinguished very quickly. At a certain 
point, fragmentation sets in, which eventually evolves into fermionization. 

It should be noted that the characteristic influence of the even atom number is on the 
transition to fermionization. Recall that this was hampered for homogeneous interactions 
owing to symmetry, which imposed that one particle had to reside near the center, x = 0. 
This no longer holds here, and in fact, Fig.lT4Tb) unveils that the spare particle is practically 
accommodated in the left well. 

The nature of the ground-state displacement is again summarized in the graph of —(x) 
(Fig. IT2|) . While the harmonic system turned out to be rather irresponsive to g , the dis- 
placement now exhibits a dramatic increase with raising go, as laid out above. It finds a 
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maximum, which corresponds to the trade-off between localizing in the left well and maxi- 
mum spreading. As before, the modulation a does not so much alter the critical g^a), but 
of course makes for a stronger maximum displacement (—x)*. The displacement decreases 
again slowly beyond that point. A notable side effect is that the displacement in the presence 
of a central barrier may in fact drop below the one without it, although of course this can 
only happen if the modulation a was smaller to begin with. That is simply because the dou- 
ble well, favoring the derealization of the atoms, not only supports the modulation's effect, 
but also accelerates fragmentation and hence — eventually — destruction of the asymmetry. 

VI. CONCLUSION AND OUTLOOK 

We have studied the numerically exact ground state of N one-dimensional bosons in 
double-well traps, where the atoms interact repulsively via a short-ranged potential whose 
strength as well as its spatial modulation could be tuned. Our approach relied on the Multi- 
Configuration Time-Dependent Hartree method, a wave-packet dynamics tool known for its 
efficiency in higher dimensions. This allowed us to study the interplay between different 
heights of the barrier separating both wells, on the one hand, and different interaction 
strengths as well as spatial modulations in all relevant regimes in a numerically exact fashion. 

For standard, homogeneous interactions, we have witnessed the transition from the weak- 
coupling mean-field regime to fragmentation and finally to a fermionized ground state for 
very large repulsion. Absent any barrier, this process requires sufficient interaction energies 
so as to compensate the one-particle (kinetic and potential) energy added by a fragmentation; 
also see [lfl| . Its signature is a broadening and eventually the appearance of N humps in the 
density profile, plus a 'correlation hole' in the two-particle density. As demonstrated, this 
also reflects in the relative occupation of the dominant natural orbital, no, which reduces 
from unity to order of l/\fN as the interaction is increased. As we turn the harmonic trap 
into a double well, then well below the barrier energy the fragmentation essentially takes 
place separately in each well, whereas way above the barrier, the situation resembles that in 
the harmonic trap. In particular, we find that for even N, fermionization is assisted, while it 
is impeded for odd N, when by symmetry one particle should be distributed over the barrier 
region. 

We have also tackled the question of inhomogeneous effective interactions, insomuch 
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as the coupling strength is assumed to be larger on one side of the trap. We have found 
the ground state to be displaced toward the side where repulsive forces are weaker. This 
displacement can be enhanced by stronger modulations, whereas the optimal interaction 
strength needed to achieve it can be lowered primarily by higher barriers. It falls off for 
stronger interactions as the fermionization destructs any imbalance effects. 

This work can be seen as a model study of the effects of inhomogeneous interactions, 
in particular a setup with an imbalance inducing a parity violation. It can in principle be 
extended to higher particle numbers, although this was not done here with an eye toward 
time and computational effort. On the other hand, there are also plenty of promising other 
configurations that come into question, such as enhancing the transmission through a barrier 
by modulating the interaction strength accordingly or using more than two wells so as to 
separate single atoms from the system. Of course, to address these problems, eventually a 
time-dependent simulation is mandatory to gain insight into realistic situations, while the 
extension to more dimensions may become inevitable. 
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